/* 
 * File:   oclFFT.cpp
 * Author: cbender
 * 
 * Created on 5. September 2014, 16:21
 */

#include "oclFFT.h"


oclFFT::oclFFT() {
    // set up everything we will use later on
    const char* cSourceFile = "CalculateFrequency.cl";
    char* cPathAndName = NULL;      // var for full paths to data, src, etc.
    char* cSourceCL = NULL;         // Buffer to hold source for compilation
    size_t szKernelLength;          // Byte size of kernel code
    gpuInput = (gpu_input_data*) malloc(sizeof(gpu_input_data));
        //Get an OpenCL platform
    ciErr1 = clGetPlatformIDs(2, &cpPlatform, NULL);
    //printf("clGetPlatformID...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clGetPlatformID\n");
        exit(0);
    }
    //Get the devices
    ciErr1 = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &cdDevice, NULL);
    //printf("clGetDeviceIDs...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clGetDeviceIDs\n");
        exit(0);
    }
    //Create the context
    cxGPUContext = clCreateContext(0, 1, &cdDevice, NULL, NULL, &ciErr1);
    //printf("clCreateContext...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clCreateContext\n");
        exit(0);
    }
    // Create a command-queue
    cqCommandQueue = clCreateCommandQueue(cxGPUContext, cdDevice, 0, &ciErr1);
    //printf("clCreateCommandQueue...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clCreateCommandQueue\n");
        exit(0);
    }
    // Allocate the Device Memory
    d_gpuInput = clCreateBuffer(cxGPUContext, CL_MEM_READ_ONLY, sizeof (gpu_input_data), NULL, &ciErr1);
    d_input = clCreateBuffer(cxGPUContext, CL_MEM_READ_ONLY, sizeof (input_data), NULL, &ciErr2);
    ciErr1 |= ciErr2;
    //printf("clCreateBuffer...\n");
    if (ciErr1 != CL_SUCCESS) {
        printf("Error in clCreateBuffer.\n");
        exit(0);
    }
    // Read the OpenCL kernel in from source file
    //printf("oclLoadProgSource (%s)...\n", cSourceFile); 
    cPathAndName = shrFindFilePath(cSourceFile, "./oclDistroFix");
    cSourceCL = oclLoadProgSource(cPathAndName, "", &szKernelLength);

    // Create the program
    cpProgram = clCreateProgramWithSource(cxGPUContext, 1, (const char **)&cSourceCL, &szKernelLength, &ciErr1);
    //printf("clCreateProgramWithSource...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clCreateProgramWithSource!!\n\n");
        exit(0);
    }
    
    //const char* options="-cl-opt-disable";
    const char* options="";
    ciErr1 = clBuildProgram(cpProgram, 0, NULL, options, NULL, NULL);
    //printf("clBuildProgram...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clBuildProgram %d!!!\n\n", ciErr1);

        size_t length;
        char buffer[32000];
        clGetProgramBuildInfo(cpProgram, cdDevice, CL_PROGRAM_BUILD_LOG, sizeof(buffer), buffer, &length);
        std::cout<<"--- Build log ---\n "<<buffer<<std::endl;

        
        exit(0);
    }
    // Create the kernel
    ckKernel = clCreateKernel(cpProgram, "calculateFrequency", &ciErr1);
    //printf("clCreateKernel (calculateFrequency)...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clCreateKernel!!!\n\n");
        exit(0);
    }
    // Create the second kernel
    //ckKernel2 = clCreateKernel(cpProgram, "fftRadix2Kernel", &ciErr1);
    ckKernel2 = clCreateKernel(cpProgram, "fftRadix2Kernel", &ciErr1);
    //printf("clCreateKernel (fftRadix2Kernel)...\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clCreateKernel!!!\n\n");
        exit(0);
    }
}

oclFFT::~oclFFT() {
    clReleaseMemObject(d_input);
    clReleaseMemObject(d_gpuInput);
    clReleaseKernel(ckKernel);
    clReleaseKernel(ckKernel2);
    clReleaseProgram(cpProgram);
    clReleaseCommandQueue(cqCommandQueue);
    clReleaseContext(cxGPUContext);
    free(gpuInput);
}

void oclFFT::setInput(input_data* input) {
    this->input = input;
}

double oclFFT::getGreatestAbsNegativeValue() {
    double result = 0;
    for (long i=0; i<this->gpuInput->numPoints; i++) {
        if(frequencyData[i].real<result) {
            result=frequencyData[i].real;
        }
    }
    return fabs(result);
}

void oclFFT::calculations() {
    __float128 carbonWeight = 0;
    __float128 integerCarbonWeight = 0;
    __float128 carbonVariance = 0;
    __float128 integerCarbonVariance = 0;
    __float128 heliumWeight = 0;
    __float128 integerHeliumWeight = 0;
    __float128 heliumVariance = 0;
    __float128 integerHeliumVariance = 0;
    __float128 nitrogenWeight = 0;
    __float128 integerNitrogenWeight = 0;
    __float128 nitrogenVariance = 0;
    __float128 integerNitrogenVariance = 0;
    __float128 oxygenWeight = 0;
    __float128 integerOxygenWeight = 0;
    __float128 oxygenVariance = 0;
    __float128 integerOxygenVariance = 0;
    __float128 phosphorusWeight = 0;
    __float128 integerPhosphorusWeight = 0;
    __float128 phosphorusVariance = 0;
    __float128 integerPhosphorusVariance = 0;
    __float128 sulfurWeight = 0;
    __float128 integerSulfurWeight = 0;
    __float128 sulfurVariance = 0;
    __float128 integerSulfurVariance = 0;
    __float128 tempMolecularVariance = 0;
    __float128 tempIntegerMolecularVariance = 0;

    // Could be done with loops but should be faster hard-coded

    if (input->numC != 0) {
        carbonWeight = (C12_MASS * input->C12 + C13_MASS * input->C13);
        integerCarbonWeight = (C12_INT_MASS * input->C12 + C13_INT_MASS * input->C13);
        gpuInput->molecularWeight += input->numC * carbonWeight;
        gpuInput->integerMolecularWeight += input->numC * integerCarbonWeight;
        gpuInput->minMolecularWeight += input->numC * C12_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numC * C12_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numC * C13_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < C13_INT_MASS - C12_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = C13_INT_MASS - C12_INT_MASS;
        carbonVariance += (C12_MASS - carbonWeight)*(C12_MASS - carbonWeight) * input->C12;
        carbonVariance += (C13_MASS - carbonWeight)*(C13_MASS - carbonWeight) * input->C13;
        integerCarbonVariance += (C12_INT_MASS - integerCarbonWeight)*(C12_INT_MASS - integerCarbonWeight) * input->C12;
        integerCarbonVariance += (C13_INT_MASS - integerCarbonWeight)*(C13_INT_MASS - integerCarbonWeight) * input->C13;
    }

    if (input->numH != 0) {
        heliumWeight = (H1_MASS * input->H1 + H2_MASS * input->H2);
        integerHeliumWeight = (H1_INT_MASS * input->H1 + H2_INT_MASS * input->H2);
        gpuInput->molecularWeight += input->numH * heliumWeight;
        gpuInput->integerMolecularWeight += input->numH * integerHeliumWeight;
        gpuInput->minMolecularWeight += input->numH * H1_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numH * H1_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numH * H2_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < H2_INT_MASS - H1_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = H2_INT_MASS - H1_INT_MASS;
        heliumVariance += (H1_MASS - heliumWeight)*(H1_MASS - heliumWeight) * input->H1;
        heliumVariance += (H2_MASS - heliumWeight)*(H2_MASS - heliumWeight) * input->H2;
        integerHeliumVariance += (H1_INT_MASS - integerHeliumWeight)*(H1_INT_MASS - integerHeliumWeight) * input->H1;
        integerHeliumVariance += (H2_INT_MASS - integerHeliumWeight)*(H2_INT_MASS - integerHeliumWeight) * input->H2;
    }

    if (input->numN != 0) {
        nitrogenWeight = (N14_MASS * input->N14 + N15_MASS * input->N15);
        integerNitrogenWeight = (N14_INT_MASS * input->N14 + N15_INT_MASS * input->N15);
        gpuInput->molecularWeight += input->numN * nitrogenWeight;
        gpuInput->integerMolecularWeight += input->numN * integerNitrogenWeight;
        gpuInput->minMolecularWeight += input->numN * N14_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numN * N14_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numN * N15_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < N15_INT_MASS - N14_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = N15_INT_MASS - N14_INT_MASS;
        nitrogenVariance += (N14_MASS - nitrogenWeight)*(N14_MASS - nitrogenWeight) * input->N14;
        nitrogenVariance += (N15_MASS - nitrogenWeight)*(N15_MASS - nitrogenWeight) * input->N15;
        integerNitrogenVariance += (N14_INT_MASS - integerNitrogenWeight)*(N14_INT_MASS - integerNitrogenWeight) * input->N14;
        integerNitrogenVariance += (N15_INT_MASS - integerNitrogenWeight)*(N15_INT_MASS - integerNitrogenWeight) * input->N15;
    }

    if (input->numO != 0) {
        oxygenWeight = (O16_MASS * input->O16 + O17_MASS * input->O17 + O18_MASS * input->O18);
        integerOxygenWeight = (O16_INT_MASS * input->O16 + O17_INT_MASS * input->O17 + O18_INT_MASS * input->O18);
        gpuInput->molecularWeight += input->numO * oxygenWeight;
        gpuInput->integerMolecularWeight += input->numO * integerOxygenWeight;
        gpuInput->minMolecularWeight += input->numO * O16_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numO * O16_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numO * O18_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < O18_INT_MASS - O16_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = O18_INT_MASS - O16_INT_MASS;
        oxygenVariance += (O16_MASS - oxygenWeight)*(O16_MASS - oxygenWeight) * input->O16;
        oxygenVariance += (O17_MASS - oxygenWeight)*(O17_MASS - oxygenWeight) * input->O17;
        oxygenVariance += (O18_MASS - oxygenWeight)*(O18_MASS - oxygenWeight) * input->O18;
        integerOxygenVariance += (O16_INT_MASS - integerOxygenWeight)*(O16_INT_MASS - integerOxygenWeight) * input->O16;
        integerOxygenVariance += (O17_INT_MASS - integerOxygenWeight)*(O17_INT_MASS - integerOxygenWeight) * input->O17;
        integerOxygenVariance += (O18_INT_MASS - integerOxygenWeight)*(O18_INT_MASS - integerOxygenWeight) * input->O18;
    }

    if (input->numP != 0) {
        phosphorusWeight = (P31_MASS * input->P31);
        integerPhosphorusWeight = (P31_INT_MASS * input->P31);
        gpuInput->molecularWeight += input->numP * phosphorusWeight;
        gpuInput->integerMolecularWeight += input->numP * integerPhosphorusWeight;
        gpuInput->minMolecularWeight += input->numP * P31_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numP * P31_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numP * P31_INT_MASS;
    }

    if (input->numS != 0) {
        sulfurWeight = (S32_MASS * input->S32 + S33_MASS * input->S33 +
                S34_MASS * input->S34 + S36_MASS * input->S36);
        integerSulfurWeight = (S32_INT_MASS * input->S32 + S33_INT_MASS * input->S33 +
                S34_INT_MASS * input->S34 + S36_INT_MASS * input->S36);
        gpuInput->molecularWeight += input->numS * sulfurWeight;
        gpuInput->integerMolecularWeight += input->numS * integerSulfurWeight;
        gpuInput->minMolecularWeight += input->numS * S32_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numS * S32_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numS * S36_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < S36_INT_MASS - S32_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = S36_INT_MASS - S32_INT_MASS;
        sulfurVariance += (S32_MASS - sulfurWeight)*(S32_MASS - sulfurWeight) * input->S32;
        sulfurVariance += (S33_MASS - sulfurWeight)*(S33_MASS - sulfurWeight) * input->S33;
        sulfurVariance += (S34_MASS - sulfurWeight)*(S34_MASS - sulfurWeight) * input->S34;
        sulfurVariance += (S36_MASS - sulfurWeight)*(S36_MASS - sulfurWeight) * input->S36;
        integerSulfurVariance += (S32_INT_MASS - integerSulfurWeight)*(S32_INT_MASS - integerSulfurWeight) * input->S32;
        integerSulfurVariance += (S33_INT_MASS - integerSulfurWeight)*(S33_INT_MASS - integerSulfurWeight) * input->S33;
        integerSulfurVariance += (S34_INT_MASS - integerSulfurWeight)*(S34_INT_MASS - integerSulfurWeight) * input->S34;
        integerSulfurVariance += (S36_INT_MASS - integerSulfurWeight)*(S36_INT_MASS - integerSulfurWeight) * input->S36;
    }

    // To compensate for the Charge, substract the weight of one electron
    gpuInput->molecularWeight -= ELECTRON_MASS;
    gpuInput->integerMolecularWeight -= ELECTRON_MASS;
    gpuInput->minMolecularWeight -= ELECTRON_MASS;
    gpuInput->integerMolecularWeightAsInt = (int) (gpuInput->integerMolecularWeight + 0.5);

    // Calculate the Molecular Variance
    // the square of the distance of the average mass of an element to each isotope,
    // summed up

    tempMolecularVariance += carbonVariance * input->numC;
    tempMolecularVariance += heliumVariance * input->numH;
    tempMolecularVariance += nitrogenVariance * input->numN;
    tempMolecularVariance += oxygenVariance * input->numO;
    tempMolecularVariance += phosphorusVariance * input->numP;
    tempMolecularVariance += sulfurVariance * input->numS;
    gpuInput->molecularVariance = (cl_real_t) tempMolecularVariance;

    tempIntegerMolecularVariance += integerCarbonVariance * input->numC;
    tempIntegerMolecularVariance += integerHeliumVariance * input->numH;
    tempIntegerMolecularVariance += integerNitrogenVariance * input->numN;
    tempIntegerMolecularVariance += integerOxygenVariance * input->numO;
    tempIntegerMolecularVariance += integerPhosphorusVariance * input->numP;
    tempIntegerMolecularVariance += integerSulfurVariance * input->numS;
    gpuInput->integerMolecularVariance = (cl_real_t) tempIntegerMolecularVariance;
}

void oclFFT::calculateMassRange() {
    // Calculate the Mass Range
    // This is the smallest Power of two that is still bigger than the
    // molecular Variance we computed, and it needs to be a power of two for the fft
    cl_ulong dPoints;
    //This is insufficient without adding the one to the end
    dPoints = (cl_ulong) (sqrt(1 + gpuInput->molecularVariance)*10)+0.5;
    gpuInput->bitSize=0;
    //printf("dPoints is %d\n", (cl_uint) dPoints);
    while(dPoints){                       // until we found the first bit.
        dPoints >>=1;
        gpuInput->bitSize++;
    }
    //printf("BitSize is %d\n", gpuInput->bitSize);
    gpuInput->massRange=1;
    gpuInput->massRange<<=gpuInput->bitSize;
    //printf("MassRange is %d\n", (cl_uint) gpuInput->massRange);
    // Enhance Resolution. Powers of 2 run well, everything else is not working
    gpuInput->bitSize += gpuInput->resolution;
    gpuInput->numPoints=1;
    gpuInput->numPoints<<=gpuInput->bitSize;
    //printf("numPoints is %d\n", (cl_uint) gpuInput->numPoints);
}

void oclFFT::reset() {
    gpuInput->molecularWeight = 0; // the average molecular weight of the molecule
    gpuInput->minMolecularWeight = 0; // the minimum molecular weight of the molecule, monoMass
    gpuInput->maxIntIsoShiftAsInt = 0; // the maximum difference in Isotope weight as Integer (can be 4,2,1 or 0)
    gpuInput->integerMolecularWeight = 0; // the average molecular weight of the molecule based on the integer values
    gpuInput->integerMolecularWeightAsInt = 0; // the average molecular weight of the molecule based on the integer values
    gpuInput->minIntegerMolecularWeightAsInt = 0; // the minimum molecular weight of the molecule based on the integer values
    gpuInput->maxIntegerMolecularWeightAsInt = 0; // the maximum molecular weight of the molecule based on the integer values
    gpuInput->molecularVariance = 0; // the molecular Variance
    gpuInput->integerMolecularVariance = 0; //the molecular Variance calculated from integer Values
    gpuInput->massRange = 0;
    gpuInput->numPoints = 0;
    gpuInput->resolution = 4;           // multiply the number of points calculated by 2^resolution
    gpuInput->bitSize = 1;
}

void oclFFT::convertMass() {
    cl_real_t mass, ratio, CorrIntMW;

    if (gpuInput->integerMolecularVariance == 0) ratio = 1;
    else ratio = sqrt(gpuInput->molecularVariance) / sqrt(gpuInput->integerMolecularVariance);
    CorrIntMW = gpuInput->integerMolecularWeight * ratio;

    for (cl_long i = 0; i < gpuInput->numPoints; i++) {
        if (i < gpuInput->numPoints / 2) {
            mass = (double) i / pow(2,gpuInput->resolution) + gpuInput->integerMolecularWeightAsInt;
        } else {
            mass = (double) (i - gpuInput->numPoints) / pow(2,gpuInput->resolution) + gpuInput->integerMolecularWeightAsInt;
        }

        mass *= ratio;
        mass += gpuInput->molecularWeight - CorrIntMW;
        frequencyData[i].imaginary = mass;
    }

}


// only copy the relevant data from the output array back into the list.
// if the weight is too big or too small for the molecule, the result is rejected
// there are also some incorrect peaks intensities from the calculation that should be zero
// and need to be rejected. if a peak is negative, it is obvious that the resulting value is nonsense.
// In order to reject the positive Values that are incorrect, the fact that there are
// as many deviations negative-wise as there are deviations positive-wise is used, and so
// the biggest negative deviation is used as a cut-off value on the positive peak values.
isotope* oclFFT::copyResult() {
    double cutoffValue = 1.1*getGreatestAbsNegativeValue(); // multiply with 1.1 to get all of them!
    isotope* result = NULL;
    isotope* temp;
    for (cl_long j = 0; j < gpuInput->numPoints; j++) {
        if ((((int) (frequencyData[j].imaginary + 0.5)) >= gpuInput->minIntegerMolecularWeightAsInt) &&
                (((int) (frequencyData[j].imaginary + 0.5)) <= gpuInput->maxIntegerMolecularWeightAsInt) &&
                (frequencyData[j].real > cutoffValue)) {
            temp = (isotope*) malloc(sizeof (isotope));
            temp->mass = frequencyData[j].imaginary;
            temp->p = frequencyData[j].real;
            temp->next = NULL;
            temp->previous = NULL;
            result = add_peak_sorted(result, temp);
            temp = NULL;
        }
    }
    return result;
}

std::chrono::nanoseconds oclFFT::calculate(isotope* out) {
    reset();
    calculations(); // calculate weights and variances
    calculateMassRange(); // calculate the mass Range
    
    frequencyData = (clcomplex*) malloc(gpuInput->numPoints * sizeof(clcomplex));
    
    // This will be done on the GPU
    // Calculate the Frequency Data

    // Allocate the Device Memory
    d_frequencyData = clCreateBuffer(cxGPUContext, CL_MEM_READ_WRITE, sizeof (clcomplex) * gpuInput->numPoints, NULL, &ciErr1);
    if (ciErr1 != CL_SUCCESS) {
        printf("Error in clCreateBuffer.\n");
        exit(0);
    }
    
    // Set the Argument values
    ciErr1 = clSetKernelArg(ckKernel, 0, sizeof(cl_mem), (void*)&d_gpuInput);
    ciErr1 |= clSetKernelArg(ckKernel, 1, sizeof(cl_mem), (void*)&d_input);
    ciErr1 |= clSetKernelArg(ckKernel, 2, sizeof(cl_mem), (void*)&d_frequencyData);
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clSetKernelArg!!!\n\n");
        exit(0);
    }
    
    // Start Core sequence... copy input data to GPU, compute, copy results back
    // Asynchronous write of data to GPU device
    ciErr1 = clEnqueueWriteBuffer(cqCommandQueue, d_gpuInput, CL_FALSE, 0, sizeof (gpu_input_data), gpuInput, 0, NULL, NULL);
    ciErr1 |= clEnqueueWriteBuffer(cqCommandQueue, d_input, CL_FALSE, 0, sizeof(input_data), input, 0, NULL, NULL);
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clEnqueueWriteBuffer!!!\n\n");
        exit(0);
    }
    
    // wait until all data is written
    clFinish(cqCommandQueue);
    // Take first time measurement
    auto start = std::chrono::high_resolution_clock::now();
    
    // Launch kernel for frequency calculation
    size_t global_item_size = gpuInput->numPoints;
    size_t local_item_size = 32;

    ciErr1 = clEnqueueNDRangeKernel(cqCommandQueue, ckKernel, 1, NULL, &global_item_size, &local_item_size, 0, NULL, NULL);
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clEnqueueNDRangeKernel!!! %d. NumPoints: %d\n\n",ciErr1, (cl_uint) gpuInput->numPoints);
        exit(0);
    }
    // wait until all data is written
    clFinish(cqCommandQueue);
    // Do the FFT
    global_item_size = gpuInput->numPoints/2;       // every Step in the FFT needs N/2 threads, each working on 2 complex data Points
    
    for(cl_int p=1; p<gpuInput->numPoints; p<<=1) { // p=1,2,4,...,numPoints/2
        //printf("Running FFT. P=%d\n",p);
        // Set the Argument values
        ciErr1 = clSetKernelArg(ckKernel2, 0, sizeof(cl_mem), (void*)&d_frequencyData);
        ciErr1 |= clSetKernelArg(ckKernel2, 1, sizeof(cl_int), (void*)&p);
        //printf("clSetKernelArg 0 - 2...\n"); 
        if (ciErr1 != CL_SUCCESS)
        {
            printf("Error in clSetKernelArg!!!\n\n");
            exit(0);
        }
        
        ciErr1 = clEnqueueNDRangeKernel(cqCommandQueue, ckKernel2, 1, NULL, &global_item_size, &local_item_size, 0, NULL, NULL);
        //printf("clEnqueueNDRangeKernel (calculateFFT)...\n"); 
        if (ciErr1 != CL_SUCCESS)
            {
                printf("Error in clEnqueueNDRangeKernel!!!\n\n");
                exit(0);
            }
        // wait until all data is written
        clFinish(cqCommandQueue);
    }
    
    // second time measurement and output of data
    auto finish = std::chrono::high_resolution_clock::now();


        
    //   Synchronous/blocking read of results
    ciErr1 = clEnqueueReadBuffer(cqCommandQueue, d_frequencyData, CL_TRUE, 0, sizeof(clcomplex) * gpuInput->numPoints, frequencyData, 0, NULL, NULL);
    //printf("clEnqueueReadBuffer (d_frequencyData->frequencyData)...\n\n"); 
    if (ciErr1 != CL_SUCCESS)
    {
        printf("Error in clEnqueueReadBuffer!!!\n\n");
        exit(0);
    }
    

    // Translate the Mass back to Real Data
    convertMass();
    
    // put the Data into the result
    result = copyResult();

    //Clean up the memory
    clReleaseMemObject(d_frequencyData);
    free(frequencyData);
    out->next=result;
    return (std::chrono::duration_cast<std::chrono::nanoseconds>(finish-start));
}